
#include "Matrix.hpp"
#include "Util.hpp"
#include "HJM.hpp"
#include "HJM_Observer.hpp"
#include "PThread.hpp"
#include "Config.hpp"
#include <stdexcept>
#include <fstream>
#include <cstdlib>
#include <unistd.h>

using namespace std;
using namespace Util;

struct Result
{
  Result() { }
  // Result(double p_sum, double p_sumsq) :
  //   sum(p_sum), sumsq(p_sumsq) { }
  int nPath;
  // short rate time series
  vector<double> sRate;
  // final test structure
  vector<double> termst;
};

struct MCArgs
{
  Random* norm;
  int nPath,nStep;
  double Tf;
  pthread_mutex_t mutex;
  HJM* iHJM;
};

// cumulative result
Result g_result;

void 
ReportResult( pthread_mutex_t mutex, Result result )
{
   pthread_mutex_lock( &mutex );

   g_result.nPath += result.nPath;
   //g_result.sum += result.sum;
   //g_result.sumsq += result.sumsq;

   pthread_mutex_unlock( &mutex );
}

void *CalcHJM( void *ptr )
{
  cout << "thread " << pthread_self() << " started" << endl;

  MCArgs* args = (MCArgs*) ptr;
  //Random* norm = args->norm;
  int nPath = args->nPath;
  int nStep = args->nStep;
  double Tf = args->Tf;
  // instance of HJM
  HJM* iHJM = args->iHJM;

  //
  // Path Generation
  //
  for ( int i=0; i < nPath; ++i )
    {
      cout << "Path " << i << endl;
      iHJM->simulate(Tf,nStep);
      Matrix<double> Rates = iHJM->result();
    }

  //ReportResult ( args->mutex, Result(nPath,sum,sumsq) );
}

int
main ( int argc, char* argv[], char* env[] )
{
  // defaults
  int nPath = 1;
  int nStep = 1000;
  double Tf = 10.0;
  int seed = 0;
  int c;
  while ((c = getopt(argc, argv, "n:p:s:t:")) != EOF)
    switch (c)
      { 
      case 'n':
        nStep = atoi ( optarg );
        break;
      case 'p':
        nPath = atoi ( optarg );
        break;
      case 's':
        seed = atoi ( optarg );
        break;
      case 't':
        Tf = atof ( optarg );
        break;
      case '?':
        cerr << "invalid option" << endl;
        exit(1);
      }

  if ( seed > 0 )
    seed = -seed;

  if ( argc - optind < 1 )
    {
      cerr << "Usage: " << argv[0] << " F0" << endl;
      return 1;
    }

  try
    {
      bool test = true;
      //
      // setup threading
      //
      PThread* PT = PThread::instance();
      int nThread = 1;
      
      Config* CFG = Config::instance();
      CFG->read("input.txt");

      // initial rates
      ifstream F0_data(argv[optind]);
      if ( F0_data.fail() )
       	{
       	  cerr << "open " << argv[optind+2] << " failed" << endl;
       	  exit(1);
       	}
      Vector<double> F0 = CSV2Vector<double>(F0_data);
      F0_data.close();

      //
      // read input data
      //
      ifstream HJM_data("data/HJM.dat",ios::binary);
      if ( HJM_data.fail() )
	{
	  cerr << "open \"data/HJM.dat\" failed" << endl;
	  exit(1);
	}
      int nX;
      HJM_data.read((char*)&nX,sizeof(int));
      cout << "Reading data for " << nX << " maturities" << endl;
      double Vtmp[nX];
      // read maturities
      HJM_data.read((char*)Vtmp,nX*sizeof(double));
      Vector<double> X(nX,Vtmp);
      // read drift
      HJM_data.read((char*)Vtmp,nX*sizeof(double));
      Vector<double> M(nX,Vtmp);
      // read eigenvectors
      Matrix<double> EigenVec(3,nX);
      for ( int i=0; i < 3; ++i )
	{
	  HJM_data.read((char*)Vtmp,nX*sizeof(double));
	  EigenVec[i] = Vector<double>(nX,Vtmp);
	}
      HJM_data.close();

      //
      //  check that everything was read correctly
      //
      ofstream ofs_vol("vol1.txt");
      if ( ofs_vol.fail() )
	{
	  cerr << "open \"vol1.txt\" failed" << endl;
	  exit(1);
	}
      ofs_vol << "Maturity";
      for ( int j = 0; j < nX; ++j )
	ofs_vol << " " << X[j];
      ofs_vol << endl;
      for ( int i=0; i < 3; ++i )
	{
	  ofs_vol << "Vol_" << i;
	  for ( int j = 0; j < nX; ++j )
	    ofs_vol << " " << EigenVec[i][j];
	  ofs_vol << endl;
	}
      ofs_vol << "M";
      for ( int j = 0; j < nX; ++j )
	ofs_vol << " " << M[j];
      ofs_vol << endl;
      ofs_vol.close();
      //
      // Initialize RNG
      //
      Normal* norm = new Normal(seed);
      pthread_mutex_t mutex = PTHREAD_MUTEX_INITIALIZER;
      // 
      // HJM simulation itself
      //

      if ( test ) {
	//
	// testing
	//
	//vector<double> XV = NRVec2Vector<double>(X);
	HJM* iHJM = new HJM(X,M,EigenVec,F0);

	// HJM observers
	vector<HJM_Observer*> Observers;
	while ( CFG->next() )
	  Observers.push_back(new HJM_Observer(iHJM,CFG->element()));

	Matrix<double> RateMean;
	// double Caplet = 0;
	// double ZCB1 = 0;
	// double ZCB5 = 0;

	cout << "\n\n*** " << nPath << " simulations ***" << endl;
	for ( int np = 0; np < nPath; ++np )
	  {
	    iHJM->simulate(Tf,nStep);
	    Matrix<double> Rates = iHJM->result();

	    // Caplet += iHJM->caplet(0.045,0.5,5.0);
	    // ZCB1 += iHJM->ZCB(1.0);
	    // ZCB5 += iHJM->ZCB(5.0);

	    if ( ! np )
	      RateMean = Rates;
	    else
	      RateMean += Rates;
	  }
	vector<double> TS = iHJM->TimeSeries();

	RateMean /= nPath;

	// Caplet /= nPath;
	// cout << "caplet = " << Caplet << endl;
	// ZCB1 /= nPath;
	// cout << "ZCB(1) = " << ZCB1 << endl;
	// ZCB5 /= nPath;
	// cout << "ZCB(5) = " << ZCB5 << endl;

	// output results
	for ( int i=0; i < Observers.size(); ++i )
	  {
	    vector<double> args = Observers[i]->args();
	    cout << Observers[i]->name() << "(";
	    cout << args[0];
	    for ( int j=1; j < args.size(); ++j )
	      cout << "," << args[j];
	    cout << ") = " << Observers[i]->value() << endl;
	  }

	ofstream ofs_test("sim.txt");
	if ( ofs_test.fail() )
	  {
	    cerr << "open \"sim.txt\" failed" << endl;
	    exit(1);
	  }
	ofs_test << "Term";
	for ( int j=0; j < nX; ++j )
	  ofs_test << " " << X[j];
	ofs_test << endl;
	for ( int i=0; i < nStep+1; ++i )
	  {
	    ofs_test << TS[i];
	    for ( int j=0; j < nX; ++j )
	      ofs_test << " " << RateMean[i][j];
	    ofs_test << endl;
	  }
	ofs_test.close();

	delete iHJM;
	exit(0);
	// end testing
      }
      // add threads
      for ( int nt = 0; nt < nThread; ++nt )
	{
	  HJM* iHJM = new HJM(X,M,EigenVec,F0);
	  MCArgs* args = new MCArgs;
	  args->norm = norm;
	  args->nPath = nPath;
	  args->nStep = nStep;
	  args->Tf = Tf;
	  args->mutex = mutex;
	  args->iHJM = iHJM;

	  // 
	  PT->AddThread ( CalcHJM, (void*)&args );
	}
      // need to gather pointers for deletion
      // collect results
      PT->join();
    }
  catch ( std::exception& ex )
    {
      cerr << ex.what() << endl;
    }
  exit (0);
}
